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ABSTRACT 


Formulas  for  curve  fitting  and  ray  computation  using 
compound  models  made  up  of  several  different  layers  of  form 
C2-Co2  (I  — )o(  Z.)  )_1  are  presented.  Examples  of  computation 

by  pocket  programmable  calculator  on  two  Sargasso  Sea  profiles, 
one  from  the  center  of  a cold  ring  eddy  are  given.  Necessary 
tables  of  the  incomplete  beta- function  and  calculator  programs 
are  Included  in  a supplement. 
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RfciY  CALCULATIONS  OF  OCEAN  SOUND  CHANNELS 
USING  A POCKET  PROGRAMMABLE  CALCULATOR 
AND  EXTENDED  FORMS  OF  THE  HIRSCH-CARTER 
MATHEMATICAL  MODEL  WITH  TABLES  OF  THE  INCOMPLETE 
BETA  FUNCTION 

L.  Baxter,  II 


I.  Introduction 

Hirsch  and  Carter1  have  given  closed  form  expressions  for  range 
and  travel  time  of  integral  numbers  of  cycles  of  ray  paths  in  the  family  of 
symmetrical  profiles  given  by: 

c1"*  <:.*“(  1 - l°*z  I **  ) (i) 

where  c is  the  speed  of  sound  at  the  vertical  distance  Z from  the  depth  at 
which  the  speed  is  and  bl  and  P are  parameters.  Pedersen  and  Gordon^, 
Weinberg^,  Stewart^,  and  others  have  developed  the  concept  of  fitting  real* 
istlc  acoustic  profiles  with  layers  of  various  curved  profile  segments  while 
matching  the  speed  and  slope  of  the  speed  at  the  layer  Interfaces.  This 
technique  prevents  the  calculation  of  "false  caustics"  and  other  artifacts 
associated  with  less  sophisticated  profile  fits  and  minimizes  the  number  of 
layers  needed  to  represent  a natural  sound  speed  profile  realistically. 

Equation  1 can  be  used  with  different  parameters  in  each  layer  of  a 
multilayer  profile  fit.  The  geometry  of  a ray  in  a layer  may  be  understood 
by  referring  to  Figure  1 in  which  a ray  from  the  reference  level  (C  * Ce  ) 
is  refracted  as  the  sound  propagates  through  higher  speed  levels,  and  Z 
(always  positive)  is  the  absolute  value  of  the  depth  difference  from  the 
reference  level.  With  Z defined  In  this  way  and  & always  positive, 


Equation  1 may  be  rewritten  as: 


fN-i 


, cT  ('-  (■*"$  ) 


The  closed  form  expressions  given  by  Hirsch  and  Carter1  for  range  and 
travel  time  apply  only  to  Integral  multiples  of  rays  from  the  reference  level 
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to  the  vertex,  , where  the  ray  becomes  horizontal.  The  portions  of 

this  path  which  we  need  to  compute  for  rays  that  traverse  several  layers 
can  however  be  expressed  almost  as  simply  In  terms  of  incomplete  beta- 
functions  which  have  been  tabulated^.  Convergent  series  for  computing 
them  directly  have  also  been  published**. 

For  profile  models  consisting  of  no  more  than  four  layers  of  this 
type  with  no  more  than  two  layers  on  each  side  of  the  sound  channel  axis, 
it  is  not  too  difficult  or  tiresome  to  do  ray  computations  with  a medium 
capacity  pocket  programmable  calculator  and  tables.  I have  done  such 
calculations,  fitting  various  natural  asymmetric  profiles  with  approxima- 
tions consisting  of  two  or  three  layers,  using  different  parameters  of 
Equation  1 in  each  layer  and  matching  speed  of  sound  and  its  derivative  at 
the  interfaces  between  layers.  In  this  paper  I outline  the  methods  and 
give  sample  results  for  two  profiles  from  the  Sargasso  Sea,  one  from  the 
center  of  a cold  eddy  and  one  outside  of  any  eddies.  I also  outline  an 
approximate  method  for  calculating  rays  that  propagate  through  a small 
horizontal  gradient  of  sound  speed  and  a method  of  calculating  range 
annoted  ray  angle  diagrams. 

II.  Incomplete  Beta-function  and  Calculator  Programs  for  Acoustic  Ray 

Computations 

Although  an  extensive  table**  of  the  incomplete  beta-function 
is  available,  the  most  important  range  of  the  variables  for  our  purpose  is 
too  sparsely  covered.  A supplement*  to  the  present  paper  tabulates  the 
necessary  detail.  The  supplement  also  contains  a Fortran  program  for 
generating  any  other  values  that  may  be  required,  and  operating  instructions 
and  listings  of  the  curve  fitting  and  ray  computation  programs  for  the 
Texas  Instruments  SR56  calculator  which  I used.  The  calculator  programs 
could  be  applied  with  little  change  to  any  equivalent  or  larger  calculator 


I 


using  algebraic  operating  system.  ^ 

III.  The  Geometry  of  Sound  Speed  Profile  Layers  in  Which  (1  - C#  ) 

We  need  the  slope  dc/dz  in  order  to  match  different  layers 

at  the  interfaces.  Differentiating  Equation  1A,  we  have: 

. \|»-l 

olt  _ at  Pc. » (<*  2)  

<*«■  X O - Z&-  <2> 

As  we  shall  show  later,  the  ray  computations  are  simpler  if  we  can 
fit  the  profile  with  layers  in  which  the  minimum  speed  of  sound  is  equal 
to  £-0  and  occurs  at  one  interface.  Therefore  the  limit  of  the  slope, 
dc/dz,  as  Z approaches  zero  is  an  important  parameter.  Remembering  that 
<*>D  and  H — & we  have  three  cases: 

Case  1.  If  y 2 O as  5T  0 , regardless  of 

the  values  of  Cc  and 

Case  2.  If  | f as  f.'+O. 

Case  3.  If  as  0 regardless  of 

the  values  of  C * and  Ot 

For  realistic  sound  speed  profiles  of  ocean  sound  channels.  Case  1 
layers  should  be  used  to  interface  at  the  axis  or  minimum  of  the  sound  speed 
profile;  the  outer  layers  may  belong  to  Case  2 or  Case  3.  This  statement 
will  be  clarified  by  the  following  more  detailed  discussion  of  layer 
geometry  for  realistic  values  of  OC  , P , C.e  and  H . 

For  refracted  rays,  the  sound  speed  does  not  usually  exceed  about 
102%  of  the  axial  speed  of  about  1.493  km/sec.  The  shape  of  the  curves  of 
Equations  1 and  2 in  the  range  of  the  parameters  for  a 2X  change  in  sound 
speed  is  most  critically  dependent  on  P , and  reasonable  changes  in  P 
may  call  for  changes  in  Qt  over  a range  of  10^®  while  changes  in  units 
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of  Z,  or  depth  variations  of  actual  profiles,  change  by  much  smaller 

ratios.  To  show  the  shape  changes  due  to  P on  the  same  axes  for  various 
values  (Figures  2 and  3) , I use  arbitrary  units  for  Z with  0**  adjusted 
to  produce  a maximum  sound  speed  change  of  about  2 X at  Z«l.  With  these 
conventions,  the  order  of  magnitude  of  is  approximately  that  which 

would  be  realistic  for  Z kilometers. 

Figure  2 shows  the  geometry  of  Equation  1 while  Figure  3 shows  that 
of  Equation  2,  i.e.  the  slope  for  the  same  values  of  the  parameters.  In 
these  figures  curves  1-6  belong  to  Case  3,  curve  7 is  Case  2 and  curves 
8-11  are  Case  1.  For  Case  1 and  Case  2 dc/dz  Increases  with  increasing  Z, 
but  in  Case  2 the  increase  is  not  significant  within  the  27.  change  in  sound 
speed.  Case  2 approximates  to  a straight  line  and  is  the  only  case  for  which 
dc/dz  at  co  is  controlled  by  the  parameter  OC  . Case  3 layers  are 
the  only  ones  in  which  dc/dz  decreases  with  Increasing  Z.  They  are 
somewhat  more  difficult  to  use  because  matching  dc/dz  to  a lower  velocity 
adjoining  layer  requires  an  interface  at  z>0.  The  process  will  be 
explained  in  the  next  section  of  this  report. 

IV.  Fitting  Ocean  Sound  Speed  Profiles  Using  Hirsch-Carter  Type  Layers 
We  can  now  see  that  the  conditions  of  Pedersen  and  Gordon^, 
matching  sound  speed  and  slope,  are  met  by  asymmetrical  profiles  (see 
Figures  4 and  5)  consisting  of  a Case  1 Hirsch-Carter  type  upper  layer 
(designated  by  U)  meeting  a Case  1 lower  layer  (designated  by  L)  at  the 
sound  channel  axis.  If  the  designations  are  used  as  subscripts  and  the 
subscript  A refers  to  the  axis  of  the  sound  channel  ca  . 

but  y -i?  tc  and  ^ ^ . The  U and  L layers  must  belong  to  Case  1 

because  only  a zero  value  of  dc/dz  at  £e  *=■  , or  minimum  sound 

velocity,  can  give  a common  tangent  at  Z-0. 
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Figur es  4 and  5 illustrate  sound  speed  profiles  from  the  Sargasso 
Sea.  The  profile  in  Figure  4 is  from  the  center  of  a cold  ring  eddy; 
that  in  Figure  5 is  from  a location  undisturbed  by  the  eddy.  To  fit 
each  of  these  with  a U and  L layer  I used  a calculator  program  which 
iterates  from  a trial  value  of  to  place  a Hirsch-Carter  type  curve 
through  C0  at  the  axis  and  points  1 ( } , ) and  2 ( ^ ~z  ^ ) 

each  marked  with  an  X.  The  dots  in  Figures  4 and  5 indicate  the  resulting 
fit. 

Where  the  fit  is  not  perfect  the  exact  values  of  OC  and  (* 
depend  of  course  on  the  chosen  points  1 and  2,  and  an  equally  good  or 
better  fit  might  appear  from  a different  choice.  It  simplifies  the  calcu- 
lations if  P corresponds  exactly  to  a tabulated  value  in  the  supplement 
or  in  Pearson^.  Therefore  it  is  worthwhile  to  try  such  a value,  as  close 
as  possible  to  the  calculated  P Ck  is  then  recalculated  to  fit 

a point  near  the  middle  of  the  layer.  The  fit  of  the  new  values  of 
and  is  checked  over  the  measured  profile.  The  values  are  adopted  if 

the  fit  seems  good  enough. 

The  fit  is  improved  if  one  does  not  try  to  cover  too  great  a range 
of  depths  with  one  layer.  The  depth  range  can  be  subdivided  into 
additional  layers  but  meeting  the  conditions  of  equality  of  sound  speed 
and  its  derivative  are  somewhat  more  complicated  when  the  interface  is  not 
at  the  sound  channel  axis.  The  sound  ray  calculations  also  become  more 
complicated  because  the  rays  must  be  divided  into  segments  that  traverse 
the  various  layers.  Inspecting  Figures  4 and  5 we  see  that  the  fit  above 
the  axis  would  be  much  improved  by  another  layer.  In  Fipures  6 and  7 
an  "M"  layer  above  the  U layer  has  been  added  to  each  of  these  profiles. 

Both  in  the  curve  fitting  and  in  the  calculations,  to  be  discussed 
later  it  will  be  useful  to  think  of  a separate  space  for  each  layer. 


In  Figure  1 the  layer  interfaces  with  adjoining  layers  may  occur  at  Q 
and  Zj  > but  the  layer  space  and  its  coordinate  system  extends 
beyond  the  portion  that  actually  fits  approximately  to  the  real 

profile.  The  ray  segments  and  Zj  Z y in  the  layer  space  are 

only  auxilliary  constructs  for  computing  the  segment  Zq  Z g which 
corresponds  to  the  real  ray  in  the  layer.  In  the  following  discussion 
the  subscripts  U,  M,  or  L are  intended  to  Indicate  the  space  in  which  a 
coordinate  is  measured. 

Figure  6 is  an  example  of  a profile  that  can  be  fitted  with  ~ I 
is  set  equal  to  at  the  chosen  interface.  at 

the  interface  is  calculated  from  Equation  2.  Then: 


<XM  = 2.  (^AAU/ 


In  Figure  7 < / . Since  OO  as  ZM  & , 

slopes  can  be  matched  only  if  the  reference  velocity,  , is  less 

than  the  velocity  C ^ at  the  interface  and  Z y\  ^ 0 at  the  same 
place. 

The  fit  was  carried  out  as  follows:  The  layer  interface  in  U was 

chosen  near  a point  of  inflection  of  the  empirical  profile.  and 

( Ac  /tfi?.)  yj  at  the  interface  were  calculated.  With  some  trial  and  error 
a layer  thickness  (2^)  max,  and  layer  parameters  P^  and 

were  selected  to  approximate  the  curvature  and  maximum  sound  speed  in  the 
empirical  profile.  (Zw)  interface  was  then  computed  from  (dc/dz)  -(dc/dz) 
interface  using  a program  based  on  Equation  2.  The  program  Iterates  from 
a trial  value  Zt  to  a more  accurate  value  of  (Zj^)  interface. 


’I 
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V.  Solutions  for  General  Ray  Segments  in  the  Hirsch-Carter  Model 


For  our  purposes  it  is  useful  to  put  the  equations  of  Hirsch 


and  Carter1  in  slightly  different  form.  Within  a layer  described  by 


Equation  1,  a ray  is  designated  by  the  angle  0O  at  velocity 
(see  Figure  1).  It  vertexes  at  2 = /T  y where 


zv  = 1 


The  variable  ^ defined  by  Hirsch  and  Carter*  as 

✓ \ I3  / ^ 


<$  = 6*z)/s<^0e 


may  also  be  expressed  at  Z by 

£ / 


5 


The  range 


, •R.x. 


covered  by  the  ray  segment  from  0 to  Z,  can 


be  written: 


I -I 


r 7» 

. 


where  B is  the  complete  beta  function  and  I is  the  ratio  of  the  incomplete 


to  complete  beta  function.  Let 


3,  =50»,  k) 


■R.,  = zv-Bt 


The  range,  ''jr  Z y.  , can  be  written 

■^hzv  = t»K0« 


,0 


The  travel  times  that  correspond  to  the  ranges  of  Equations  8 and  9 may 
be  written 


T - T* 

I n?  — 


T«, 


c~i-  f~ j ~^) 

cnnr  (I  - >“  a.  5^7-  j 

= — { l - 7 

e.  c«£„  L (KO*.  5 


The  values  of  the  complete  beta  function,  and  B2  are  constants 
for  a given  layer  of  a profile.  They  may  be  calculated  or  taken  at  once 
from  the  tables.  The  relative  values,  1^  and  I2,  of  the  Incomplete  beta 
function  depend  on  0o  and  Z through  Equations  4 and  6,  and  the  tables. 

The  range  and  travel  time  of  any  segment  such  as  Q-S  (Figure  1)  is 

easily  obtained  as  a difference  between  values  computed  by  the  above  equations. 

Programs  for  Equations  4,  6,  8,  9,  10  and  11  fit  easily  in  the  SR56 
calculator  when  1^  and  I2  are  entered  from  tables.  Note  that  if  H - 2 y,  , 
Ij  and  I2=l  and  Equations  8 and  10  reduce  to  those  given  by  Hlrach  and 
Carter1.  To  obtain  total  ranges  for  N axis  crossings  the  values  can  of 


course  be  multiplied  by  2N  as  is  done  by  Hirsch  and  Carter. 


-9- 


I I 


VI.  Calculation  va  Axial  Angle  of  Range  and  Travel  Tine  at  the  End  of  Loops 
Above  and  Below  the  Sound  Channel  Axis  and  at  the  End  of  a Complete 
Cycle 

Since  the  classical  ray  acoustics  paper  of  Ewing  and  Worzell^, 
sound  channel  computations  have  often  been  presented  by  plots  of  range 
and  travel  time  of  loops  above  and  below  the  axis  and  of  a full  ray  cycle, 
all  vs  axial  angle  as  the  independent  variable.  These  data  are  presented 
for  the  three-layer  fits  to  the  eddy  and  Sargasso  Sea  profiles  in 
Figures  8,  9,  10  and  11.  The  procedure  for  calculating  these  plots  is 
described:  first  for  the  simpler  case  of  Figures  8 and  9 where  s / , 
and  then  the  modifications  for  Figures  10  and  11  where  ^ /. 

Axis  to  axis  loops  that  do  not  penetrate  into  a second  layer  are 
computed  by  straight-forward  application  of  Equations  4,  8 and  10  with 
the  factor  2N  equal  to  2,  Ij  and  I 2 equal  to  1,  and  @e  equal  to  the 
axial  angle,  . On  the  lower  side  of  the  axis  where  the  profile 

fit  has  no  second  layer,  the  full  range  of  axial  angles  may  be  covered  this 
simply. 

The  Z coordinate  of  the  interface  in  U layer  space  can  be  written 
. When  y')  becomes  greater  than  i~£.  ^ , the 


calculations  can  be  simplified  if  is  chosen  so  that  ^ is 

exactly  equal  to  a value  of  X printed  in  the  tables  of  the  incomplete 


beta-function.  Omitting  the  subscript  U: 


s Zt  /x 


(12) 


The  axial  angle,  C7^  , for  this  ray  is  given  by: 


= *c  *8*[(a*v)<’A]  (U) 

When  | , the  reference  sound  speed,  , for  the  M layer  is 


equal  to  , the  sound  speed  at  the  Interface.  The  reference  angle 


in  the  M layer  Is  calculated  by 


C cs 


With  and 


tabulated,  one  returns  to  the  program  for 


Equations  4,  8 and  10.  Taking  1^  and  I2  directly  without  interpolation 
from  the  tables,  range  and  travel  time  for  the  portion  of  the  loop  in  the 
U layer  is  computed  using  the  Q just  found.  Ij  and  l£  in  the  layer 
equal  1 in  this  case.  The  portion  of  the  ray  in  the  M layer  is  computed 
using  the  (B.)  M equivalent  to  Q ^ from  Equation  14.  The  values  of 
range,  travel  time,  and  distance,  Z-y  , in  the  U and  M layers  are  added 
to  obtain  the  values  plotted  for  the  ray  loop  above  the  axis.  These  range 
and  travel  time  values  are  added  to  those  computed  for  the  seme 
below  the  axis  to  obtain  the  values  for  the  full  ray  cycle. 

When  < | , as  in  Figures  10  and  11,  the  segments  in  M must  be 
computed  differently.  and  Equations  9 and  11  must  be 

used  instead  of  8 and  10.  Ij  and  I2  in  the  M layer  are  not  equal  to  1. 

One  could  use  directly  the  that  correspond  to  &K  by 

Equation  14,  but  one  would  have  to  interpolate  in  the  tables  for  1^  and  Ij. 

It  is  easier  to  defer  the  interpolation,  doing  it  in  r ^ 

a later  stage.  This  is  done  by  using  Equations  12  and  13  on  the  M layer 
after  they  have  been  used  on  the  U layer.  Equation  13  however  is  understood 


IQ)  = SVC  SI* 


MMMH MMH 
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where 


Is  the  value  of 


/|^l^  xd  buv  vuiuu  wa  ^ ' 4;  thit  correspond!  to  e 

tabulated  value  of  X and  not  to  . Range,  travel  time,  and  ^ ^ 

computed  in  the  M layer  for  ( % are  interpolated  to  find  the 
values  for  ( ^ that  do  correspond  to  @ ^ . The  results  are  added 

to  those  for  the  U layer  as  before. 

VII.  Calculation  of  Arrival  Times  for  the  Eigen  Rays  for  a Source  and 
Receiver 

This  problem  is  merely  an  extension  of  the  techniques  used  in 
Section  VI.  First  one  adds  the  appropriate  segments  to  obtain  a plot  of 
range  vs  axial  angle  for  the  source  and  receiver  depths  and  the  possible 
types  of  path.  Figure  12  is  an  example  of  this  step.  One  interpolates  to 
find  the  axial  angles  of  each  path  at  a given  range.  Then  range  and  travel 
time  are  computed  for  these  axial  angles.  Due  to  limited  precision  in  the 
first  interpolation  the  ranges  will  differ  slightly  but  the  average  sound 
speeds  will  be  correct  for  each  path  at  the  desired  range.  A second 
linear  interpolation  will  adjust  all  the  travel  times  to  the  correct  range. 
Table  I illustrates  the  result  at  a range  of  705  km  in  Figure  12  and  rays 
of  order  14  through  16. 
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TABLE  1 


Travel  Time  at  705  km  of  rays  of  order  14,  15  and  16  In 
Sargasso  Sea  Profile 


AT 


N 


7.16964 

471.65950 

0.0000000 

16 

7.21127 

471.65254 

0.0069642 

16 

7.40994 

471.604778 

0.054726 

16 

7.46963 

471.59383 

0.1001213 

16 

8.29557 

471.438324 

0.2211799 

15 

8.40314 

471.414127 

0.2453772 

15 

8.57449 

471.36485 

0.2946512 

15 

8.69705 

471.33525 

0.324255 

15 

9.70808 

471.11683 

0.54267 

14 

9.87478 

471.068211 

0.59129 

14 

10.03846 

471.015474 

0.64403 

14 

10.21703 

470.959929 

0.69958 

14 
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VIII.  Calculation  of  the  Relative  Intensity  or  Focusing  Factor 
a.  Relative  Intensity  Except  at  Caustics 

Q 

Brekhovskikh  defines  a "focusing  factor"  f»l/I0,  the  ratio 
of  the  acoustic  intensity  I at  a given  point  in  the  homogeneous  medium  to 
the  acoustic  intensity  IQ  in  a homogeneous  medium  at  the  same  distance. 

He  shows  that  when  H and  the  point  is  not  a caustic. 


f = 'R /sin 


(16) 


ft 


where  is  the  horizontal  angle  at  the  given  point  and  the  derivative 

is  evaluated  for  the  ray  that  passes  through  the  point. 


P?  ~ J\TC  CPS 


(*  tc.') 


(17) 


) may  be  obtained  graphically  as  the  slope  from  a 
plot  like  Figure  12. 

b.  Relative  Intensity  at  a Caustic 

in  Figure  12  the  four  rays  of  a given  order  appear  in  two  pairs. 
Each  pair  appears  to  join  at  a point  for  an  axial  angle  slightly  less  than 
7°.  The  scale  is  too  coarse  to  show  the  detail  in  the  neighborhood  of  the 
supposed  point  which  is  really  the  location  of  a caustic.  Figure  13  shows 
the  "point"  of  the  lower  pair  of  order  17  on  a greatly  expanded  scale.  The 
method  of  calculating  this  detail  will  be  discussed  after  I outline  its 
application. 

We  have  been  interested  in  comparing  the  relative  intensity  of 
caustics  in  differing  profiles  but  at  a given  range.  Although  ordinary  ray 
theory  fails  at  a caustic,  Brekhovskikh^  discusses  a method  of  calculating 


intensity  at  a caustic  from  ray  parameters.  The  full  expression  involves 
an  Airy  function  and  is  rather  complicated,  but  to  compare  the  maxima  of 
caustics  under  different  conditions  without  computing  the  true  relative 
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); 

p 


intensity  at  any  point  the  expressions  can  be  shortened.  In  the  notation 
o£  this  paper,  and  discarding  factors  that  don't  vary  much  in  actual  sound 
velocity  profiles,  relative  intensity  at  a given  large  range  and  a given  ^ 

acoustic  frequency  is  inversely  proportional  to  tan  0^  sin  ^ 

£)  r 
where  (see  Equation  17)  is  the  angle 

with  the  horizontal  of  a ray  tangent  to  the  caustic.  The  method  of  computing 

the  data  for  Figure  13  enables  us  to  evaluate  the  derivative 

the  other  factors  are  obvious. 

To  calculate  the  range  of  a ray  near  the  caustic,  we  measure,  in 
Figure  10,  the  slope,  Sf,  and  intercept  Ip,  of  the  full  cycle  ray  that 
vertexes  at  the  receiver  depth.  The  values  are:  Sp*3  km/degree,  If=22  km. 

We  measure  also  the  slope,  Su,  and  intercept  Iu  of  the  upper  branch. 

Su=0  km/degree.  Iu=10  km.  We  then  calculate  range  from  the  vertex  vs 
axial  angle  for  segments  to  the  receiver  depth  of  rays  that  vertex  slightly 
shallower.  The  result  appears  as  Figure  14.  We  note  that  the  range 
increment,  rx,  due  to  this  segment  is  approximately  the  parabola 


(18) 


where,  in  the  given  example,  K-2.9781  and  ( & , the  axial  angle  of 
the  ray  that  vertexes  at  the  receiver  depth,  *6.892  degrees. 

Let  the  angular  difference^  , - MrJ  = r . Total 


range 


of  a ray  of  order  N in  the  vicinity  of  the  vertex  can  be  written  as  follows: 


n - q(t*+ V'fK  y)^  k (f 


(19) 


where  Q = 3/2  or  1/2  depending  on  whether  there  is  or  is  not  an  extra 
upper  loop  in  the  group  of  rays  under  consideration.  Figure  13  is  a plot 


■5M 


of  Equation  19.  As  Indicated  by  Brekhovskikh^  the  cauatic  occura  where 
— 0 , on  the  branch  of  the  curve  with  the  minus  algn.  Now: 

oTR/$^  = - Qi Su  — K '///Z  Cf  ^ (2 


Therefore,  at  the  caustic 


cf  = K/4-  (Q‘Su  + N-Sf)  (21) 

Differentiating  Equation  20,  we  have 

= K x cf  + 

' A \ S (22) 

Evaluating  Equation  22  at  the  cauatic  by  substitution  of  Equation  21,  we 
find  that  ^ 

«*-  (G'Su  + 

V N (23) 

IX.  Calculation  of  New  Axial  Angles  of  a Ray  that  Propagates  from  one 
Profile  to  Another 

Milder10  has  shown  that  if  the  change  from  one  profile  to 
another  is  sufficiently  gradual,  there  is  an  invariant  called  the  character- 
istic time.  This  invariant  can  be  calculated  by  the  equation: 

r=  (t-x.cos  ^/ca  )/Zt 

V (24) 

where  J is  the  characteristic  time,  T is  the  full  cycle  travel  time,  X is 
the  full  cycle  range,  ^ is  the  axial  angle  of  the  ray  and  f\  is 
the  speed  of  sound  at  the  axis.  The  conditions  on  the  horisontal  gradient 
for  validity  are  given  in  detail  by  Milder  for  both  wave  and  ray  theory. 

In  ocean  sound  channels  a horizontal  gradient  as  small  as 
.03  m/sec/km  is  safe. 


Figure  15  is  a plot  of  characteristic  time  vs  axial  angle  for  the 


profiles  that  we  have  been  considering.  To  find  the  angle  in  one  profile 
that  is  equivalent  to  an  angle  in  another  one  finds  the  value  J corres- 
ponding to  the  angle,  e\  , in  the  first  profile  moves  horizontally 
to  the  curve  for  the  second  profile  and  under  the  same  J one  finds  the 
value,  <9/^  , in  the  second  profile. 

X.  Calculation  of  Range  Annotated  Ray  Angle  Diagrams 

Flatted  and  Cox*^  have  discussed  the  range  annotated  ray  angle 
diagram  and  its  applications.  A program  (for  the  pocket  programmable 
calculator)  adapted  to  generating  data  for  such  a diagram  is  included  in 
the  supplement.  The  depth  difference  Z from  the  reference  level  and  the 
angle  & are  computed  from  the  range  on  a segment  shorter  than  that 
from  reference  level  to  vertex.  Longer  paths  are  plotted  by  symmetry 
and  addition.  There  are  two  cases:  One  where  the  given  range  is  R„ 
in  Equation  8,  the  other  where  the  given  range  is  in  Equation  7. 

In  the  first  case: 


(25) 


r,  = i\z  PCt^ej/y,  z: 


In  the  second  case: 


r = i -k 


2i. 


B,Z, 


(26) 


The  value  of  1^  is  used  to  obtain  X from  the  Ix  tables  of  the  incomplete 
beta -function.  Then 


Z.=ZV* 


(27) 


and  from  Equation  1 


t/cm  = 0 - (etEl  ) 


(28) 


0 
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Q - 4*c.  crs  (t  c-os  9*/c*) 

0 - A*'  e.j  («-«  &c/ ) 


XI.  Notes  on  the  Values  of  P In  Asymmetric  Profiles  Based  on  the 


Hirsch-Carter  Model 


Hirsch  and  Carter  have  pointed  out  that,  in  symmetric  models 


of  the  near  axis  sound  transmission,  the  observed  time  dispersion  of 
arrivals  occurs  only  in  that  subset  of  the  P family  for  which  I ^ ^ 


The  actual  sound  channel,  however,  is  grossly  asynmetrical.  Because  the 


refraction  below  the  axis  is  so  much  weaker  than  that  above,  rays  at 


more  than  a very  small  axial  angle  will  spend  much  more  time  below  the 


axis  than  above  it  so  that  the  overall  dispersion  pattern  is  like  that  of 
a symmetric  channel  with  P near  the  below  axis  value  of  1.25  or  1.26, 
the  profiles  above  the  axis  are  fitted,  however,  with  values  of  P 


between  2 and  3.  This  would  tend  to  reduce  the  dispersion  below  what  one 


would  get  by  reflection  of  the  lower  half  of  the  channel. 
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' PTIONS  TO  FIGUBES 


Figure  1 Path  of  a i from  reference  level  to  vertex  within  a layer. 

Calculation  of  range  and  travel  time  of  the  segment  from 
0 to  Z or  that  from  Z to  Zv  is  discussed  in  the  text.  The 
general  segment  Q to  S can  be  expressed  as  a difference  of 
segments  of  either  of  the  above  types. 

Figure  2 Geometry  of  sound  speed  profiles  described  by  Equation  1A. 

C0=l. 49275  in  all  the  curves.  The  other  parameters  are 
listed  after  the  indicated  number  of  each  curve  as  follows: 

-t« 

1.  <*  * 10  j ; 

2.  * • 3 -*0  ; 

3.  oc.  = 3 X »0“4  J far  2*5  o 

4.  <x  * x {*  • 

5 r l.tX  ,0-\  (3  s .<50  J 

6 <*  = 7>o * , r * ^ i 

7.  04=  .O'*  , (**  1.0 

8.  .m  , f8,1'5  > 

9.  ,1Q  , ; 

10.  < * ‘MS;  *****  * 

a.  * - •*'<*  , <'=  3'°  • 

Figure  3 Slopes  of  the  sound  speed  profiles  of  Figure  2 on  a logarithmic 
plot. 

Figure  4 A sound  speed  profile  from  the  center  of  a cold  ring  eddy  is 
indicated  by  the  continuous  line.  Two  layers  according  to 
Equation  1 have  been  fitted  at  the  axis  and  the  points  marked 
by  X.  The  calculated  sound  speeds  at  other  points  are 
indicated  by  dots.  The  parameters  are:  ©<  w x ft  — 2. *^7*1 4U 

Ce—  I • *1  ^ *4  fT)  / d , O *3^(7  % j 1mc| 

Figure  5 A sound  speed  profile  in  the  Sargasso  Sea  outside  of  the  eddy 
is  indicated  by  the  continuous  line.  Two  layers  according 
to  Equation  1 have  been  fitted  at  the  axis  and  at  the  points 
marked  by  X.  The  calculated  sound  speeds  at  other  points 
are  indicated  by  dots.  The  parameters  are:  .* 

(J, -=  x,t><rse>e. 


i 
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Figure  6 The  upper  portion  of  the  eddy  profile.  Figure  4,  is  fitted 
by  two  layers.  The  previous  fit  is  retained  below  the 
sound  channel  axis.  The  new  parameters  are  as  follows: 

«XW  ® ^ , (cJ)v  ~ f'***S£7 ; = 0.  CC| 

5 I .COCO  ^ | 00 

Figure  7 The  upper  portion  of  the  Sargasso  Sea  profile.  Figure  5, 

is  fitted  by  two  layers.  The  previous  fit  is  retained  below 
the  sound  channel  axis.  The  new  parameters  are  as  follows: 

» 0.31  S',  pu  r 1,000,  (0„  * l.tltlf,  «/<TjS 

*n*c.>i>0 -(z.\  ,(•«* ,.««  =o.oiife 

Figure  8 Range  and  vertex  depth  vs  axial  angle  for  the  three  layer 
fit  of  the  eddy  profile  (Figures  4,  6).  Range  of  a loop 
above  the  axis,  one  below  the  axis,  and  a full  ray  cycle  are 
shown.  Vertex  depth  is  for  the  upper  loop  and  is  therefore 
the  shallowest  point  reached  by  the  ray. 

Figure  9 Travel  time  vs  axial  angle  for  the  same  profile,  fit,  and 
paths  as  Figure  8. 

Figure  10  Range  and  vertex  depth  vs  axial  angle  for  the  three  layer 
fit  of  the  Sargasso  Sea  profile  (Figures  5 and  7) . The 
data  presented  corresponds  to  that  presented  in  Figure  8. 

The  constancy  of  range  of  the  upper  loop  over  the  axial 
angles  0-10.4  is  a property  of  the  fit  with  = 2.0  as  noted 

in  the  Hirsch-Carter*  paper. 

Figure  11  Travel  time  vs  axial  angle  for  the  same  profile,  fit,  and 
paths  as  Figure  10. 

Figure  12  Range  vs  axial  angle  of  high  order  rays  in  the  Sargasso  Sea 
profile  (Figures  5 and  7).  The  order  is,  of  course,  the 
number  of  loops  below  the  axis.  The  receiver  depth  is  .85  km. 

The  source  is  on  the  axis.  There  are  four  rays  belonging  to 
each  order.  This  figure  may  be  used  as  described  in  the  text 
to  find  axial  angles  of  eigen  rays  at  a given  range. 

Figure  13  Range  vs  axial  angle  for  two  rays  of  the  17th  order.  This 
figure  demonstrates  the  formation  of  a caustic  as  discussed 
in  the  text. 


Figure  14 


Figure  15 
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Range  vs  axial  angle  of  a ray  segment  from  vertex  to 
receiver  depth  (.85  km,  i.e.  .4  km  above  the  axis)  in 
the  Sargasso  Sea  profile  (Figures  5 and  7).  The  solid  line 
shows  values  calculated  using  the  Hirsch-Carter  model  with 
tables  of  the  incomplete  beta-function.  The  circles  show 
values  from  the  parabolic  fit,  ^ C } p 

n hi  — ^ 


with  K r 2 A 1 9 \ 


and  = b.W-L' 


Characteristic  time  vs  axial  angle.  This  figure  can  be  used 
for  estimating  the  change  in  axial  angle  when  a ray  propagates 
through  a transition  region  from  one  sound  velocity  profile 
to  another. 
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Supplement  to 


RAY  CALCULATIONS  OF  OCEAN  SOUK)  CHANNELS 
USING  A POCKET  PROGRAMMABLE  CALCULATOR 
AND  EXTEIDED  FORMS  OF  TIE  HIRSCH-CARTER 
MATHEMATICAL  MODEL  WITH  TABLES 
OF  THE  INCOMPLETE  BETA-FUNCTION 


Tables  of  the  Incomplete  Beta-function 


The  complete  beta-function 


W r(?K(^/r  (p+'O 


is  given  by 


The  incomplete  beta-function 


^'sometimes  called  the  relative  incomplete  beta 


The  function 


function  is  given  by 


In  these  tables  and  in  those  of  Pearson  is  given  as  a function 

of  and  y.  , while  is  tabulated  at  the  top  of  each 

column  of  JL*  Cv  «v  In  the  present  tables  CJ  is  taken  equal  to  .5, 

> * 

the  only  value  it  assumes  in  the  Hirsch-Carter  model  equations. 

In  the  paper  to  which  this  supplement  is  appended  the  Hirsch-Carter 


model  equation  is  written 


and  it  is  shown  that  range  and  travel  time  can  be  written  for  sound  ray 
segments  in  terms  of  the  following  quantities 


•z-  s 


3,  = ,'5) 

2,  ” Z* 

3*.=  3(>*,.*'> 


where 


(6) 

(7) 

(8) 

(9) 

(10) 


and  ^ v is  t*ie  value  of  ~Z—  for  which  a sound  ray  vertexes  (i.e. 
becomes  horizontal  at  its  maximum  or  minimum  depth  of  excursion) . To  use 
these  tables  for  sound  ray  calculations,  3 ( is  taken  from  the  top  of 
the  column  with  -p  — ; X ( is  taken  opposite  X from  the 

same  column;  3 j_  is  taken  from  the  top  of  the  column  with  p - / *+■ 
and  T-j_  is  taken  opposite  ^ from  that  column. 

For  some  computations  it  may  be  necessary  to  enter  the  tables  with 
a value  of  X.  and  interpolate  to  find  a value  of  0( 
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B.  Fortran  Prograns  for  Calculating  Tables  of  the  Incomplete 


Beta-Function 


Abramowitz  and  Stegun2  give  a series  expansion  for  the  function 

which  is  equivalent  to  the  following: 

T / V-  «Wfi  i r 0 oc" 

* ^ f J_  <u> 

This  converges  well  if  ^ ^ . 

For  • 5”  ^ I the  symmetry  relation: 


may  be  used  to  evaluate  Equation  11  within  its  region  of  good  convergence. 
The  following  Fortran  programs  use  equations  1,11,  and  12  and  a polynomial 
approximation2  to  the  garana  function  to  tabulate  T*  • 

The  main  program  "TABLE"  calls  the  other  functions,  accepts  the  input 
parameters,  and  prints  the  output.  When  TABLE  is  called  it  requests  a 
logical  unit  number  for  the  output  with  "ENTER  IPRNT".  When  this  is  typed 
in,  the  line  "ENTER  A,  B,  N,  XO,  K"appears.  Type  in,  in  free  field  form, 
the  quantities  p,  q,  N,  XO,  K where  p and  q are  the  quantities  we  have  been 
using  (p  and  .5),  N is  the  number  of  terms  in  the  expansion  (25  for  6 place 
accuracy) , XO  is  just  smaller  than  the  smallest  desired  value  (X0=0  for  a 
complete  table  column) , and  K equals  the  number  of  tabular  entries  desired 
in  the  column  to  be  printed.  When  these  are  typed  in,  the  line  "ENTER  DEUI" 
appears.  Type  in  the  increment  between  successive  values  of  "X  and  the 
program  will  type  out  the  table  with 

* - XO  + X *3>ELX 

with  X — I ft?  K. 

Finally,  the  line  "NEW  VALUES,  YE  or  NO"  appears.  Type  in  NO  to  exit  the 
program  or  YE  to  repeat. 


FTN4*  L 

PROGRAM  TABLE 

5005  FORMATC "ENTER  IPRNT")  

5000  FORMAT  < “ENTER  A*  B*  N*  XO*K“> 

-5012  F0RMAT<5X/I4, 10X,F14.'7>-  ' ' 

5010  F0RMHTC5X,  FA.C..H-X, 

— 5020~  FORMAT  < 7XT”  X "7  20X/  M I M > "l: — ’ 

5020  FORMAT < "NEW  VALUES*  VE  OR  NO") 

— - 5040  FORMAT  <A2)  *.  '.T~T 

5045  F0RMAT<7X* "A"* 20X* "B">  ' v\ 

5047“  FORMAT  C7X7  * N "*~20X*‘  “ K"-)  ' —T-rT 

5060  FORMAT<  "ENTER  DELX”)  :X  K 

IVES=2HYE~ 

IITTV-1  ' ".  •■■■  . 

— “1 OTT  V=1  — ' — “ 

MRITE<IOTTV*5005)  /<•  . 

—  READCT I TTV7  V )T  PRNT  ' *~“ 

100  CONTINUE  , ?. 

“ UR ITECI OTT V7 5060)  " 

READ<IITTV*  *>DELX 

—  MR  ITE’CT  OTT  Vi  5000)  

READ<IITTV,  *>A,B,N*  XO,K 

HRTTEC I PRNT7 5045 ) " ~ 

MRITECIPRNT* 5010>A* B * 

—  -WRITEaPRNT;  5047)  T" ':V..  .•••/ 

MRITE<IPRNT*  5043>N,  K > ^ 
”"5048  -pORMATT5Xn4710X*  I14>' 

!->  .■  MRITECIPRNT* 5020)  .viA 

— T — : — BO~"200— I**1*TC — ■;  -V . ti-u-  ' 

v:  '-  ■ X»XO+DELX*I  '’v!  . 

— — cotsbi  <x;  a;  B7W>  : - 

; ' MRITECIPRNT*  5010)X*  OUT ' T;!' 

"'“‘T200  CONTINUE ~ 

HRITE<IOTTV*  5030) 

— READ  Cl  I TTV7 5040 FI QUER  ~ T’ 

IF<IQUER.  EQ.  IVES)GO  TO  100  > 


FUNCTION  BI<X,A,B,N) 

X*=X**A 
Xb=<1.0-X)**o 
FCrh=XH*XB/(«*bETA(A,B)> 
IF(X.GT.0.50)FCIh=XA*XB/(B*BElA(A,B) ) 
aiN=r;0  - - - 

AP=H+1.0 

~IF(X.GT70»5O)APsB+l .0  

Adsm+o 

— x**=x — ---■ — 

IF(X.GT.0.50>Xu=l.0-X 

" DO  10  1=1, N ” - 

/ XI=XQ**I 

/ Yi=r 

' BE 1=BETA  CAP , YI ) 

BE2=BETA(AB, YI ) 

BI N=B  I N+X I *&E 1 /BE2 

II'  CONTINUE  - “ — : 

BI=BIN*FCTR  • - ■ 

IF(X7GT70.50)B1=  1 .0-BI — * ' 

RETURN 

’ " END : 


FUNCTION  feETA(P,Q) 

'PG=GAMMA<P)  ~ 

QG=GANNA(U)  ■ . . 

-pQG=GANNACP+Q> ' 

BETA=PG*QG/PuG 
“RETURN  "~T. _ ‘“ 

END  . ; 

~ 

PUNCH  ON  GANNA  (A)  { 

- DOUBLE  PRECISION  AKC8) ,RG, AD, ALD,DJ,AI, AG  ■ ! 
DATrt  AK/-.577191652D0, .9SS205891D0, 

1 -. 89705695 7D0, .9 18206857 D0 ,-.75670407800 , 

'—2  .482199394D0,- . 1 955278 18D0, .05586854500/ 

AD=a 

■ ag=i.0  : ••• 

KG=1.0 

”~if<ad;ge.'K0D0)go  to  19  " " ' 

do  10  1=1,8 

: " RG=8Gl-AD**I*AKCI>  T 

10  CONTINUE 

rg=r  g/ad‘  ; : : * 

CONTINUE 

“15  GAKMA=RG~"7  "T"  •••-•.  - 

RETURN  ’•  ' ■’  V. 

19  j=ad  : — 7 : — . — — — - 

SW>\  IFCAD.EQ.1.0D 0>J=0  V/!-' 

■ I F <A  D~.Efcr.T2;  0D0 ) J=  1“  " 

' DJ=J  .•  Vij-:  -': 

- ALD=aD=D3" 1 7“ ” 

; — DO  25  1=1,8  ' ; >■ 

’ 7”  ~RG=R  GVA  LD**I *A  K C I ) ~ ' * - ‘ i 

“25 "CONTI NUE“~  ~ • - - — 

ifcad.le.2.0D0)go  to  15  

^26  CONTINUE  r-v-  •' V:  • _ . 

a g=ag»cald+ad  ; - ~:;1  • ___  _ _ _ / 

a‘  * IFCJ.GT.0JGO  TO  26  •-  ...  7 

“~fcG=AG*AG“;'  *~r;  ; 

7 GO  TO  15 

— end'  “ — ——7—— — — - 

• ends  ' •••  • 


r 
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C.  Programs  for  Profile  Fitting  and  Acoustic  Ray  Plotting  in  Long 
Range  Deep  Ocean  Sound  Transmission  Studies 

The  following  programs  written  for  the  Texas  Instruments  SR56 
programmable  pocket  calculator  would  be  easily  adaptable  to  other 
programmables  using  algebraic  notation. 

Tables  CIA  and  C1B  give  operating  instructions  and  a listing  of  a 
program  that  finds  values  of  & and  P for  a Hirsch-Carter  model  fit 
to  the  three  points:  Cv^  O ; and  "2.  . The  program 

iterates  from  a trial  value  of  P to  find  OC  and  d more 

accurate  value  of  P . The  iteration  can  be  continued  to  any  desired 
precision.  The  operating  steps  7 and  8 use  the  final  values  of  and  P 

calculate  C as  a function  of  Z.  for  comparison  with  the  empirical 
profile. 

Tables  C2A  and  C2B  give  operating  instructions  and  a listing  of  a 
program  that  calculates  slope  ( and  sound  speed  as  a function  of 

Z in  a Hirsch-Carter  type  profile  with  parameters  , P , and  C*  t . 

Tables  C3A  and  C3B  give  operating  instructions  and  a listing  of  a 
program  that  calculates  the  depth  increment  Z that  corresponds  to  a given 
sound  velocity  in  a Hirsch-Carter  type  profile  with  parameters  , p* 

and  C 

Tables  C4A  and  C4B  give  operating  instructions  and  a listing  of  a 
program  that  calculates  the  depth  increment,  Z,  that  corresponds  to  a given 
slope  dc/dz,  in  a Hirsch-Carter  type  profile  with  parameters  , P'  , 

and  £*o  . The  program  iterates  from  a trial  value  of  Z to  find  a more 
accurate  value.  The  iteration  can  be  continued  to  any  desired  precision. 
Tables  C5A  and  C5B  give  operating  instructions  and  a listing  of  a 


to 


program  that  computes  range  and  travel  time  of  a ray  segment  or  multiple 
ray  segments  in  a Hirsch-Carter  type  profile.  See  Section  V of  the  basic 
paper  for  more  detail  of  the  equations  programmed.  The  parameters 
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which  for  a case  1 profile  (see  Section  IV  of  the  basic  paper)  fitted  at 
the  axis  is  the  axial  angle  . Range  may  be  specified  from  an  axis 

crossing  or  a vertex.  Tables  of  the  incomplete  beta-function  are  used 
to  convert  X(  ~ X*  & V)  to  0(  . This  program  can  be  used  to 
generate  data  for  a range  annotated  ray  angle  diagram  as  described  by 
Flatte^  and  Cox^. 
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-^Formulas  for  curve  fitting  and  ray  competition  using  compound  models 
made  up  of  several  different  layers  of  form  ( f ) f * are 

presented.  Examples  of  computation  by  pocket  programmable  calculator  on 
two  Sargasso  Sea  profiles,  one  from  the  center  of  a cold  ring  eddy  are  given. 
Necessary  tables  of  the  incomplete  beta-function  and  calculator  programs 

are  included  in  a supplement. 
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